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A semiclassical approximation to the time evolution of coherent states may be derived from a 
saddle point approximation to the exact quantum propagator, and in general it involves complex 
classical dynamics. We generalize previous one-dimensional results to d dimensions, and for the case 
d = 2 we present several applications. We also consider other simple approximations that depend 
only on real classical trajectories, but are not initial value representations. These approximations 
are able to reproduce interference and tunnelling effects, and involve propagating a few classical 
initial conditions compatible with the quantum uncertainties. 



I. INTRODUCTION 



Semiclassical propagators involving complex classical 
trajectories in real time have appeared in the coherent 
states representation around 25 years ago |3j ■ A sta- 
tionary phase approximation to the transition amplitude 
(z'|e~*^"^/''|z), where \z) is a coherent state, leads to tra- 
jectories satisfying the usual Hamilton equations subject 
to special boundary conditions that can only be satisfied 
in a complexified phase space. Numerical calculations in 
this representation have been done both for chaotic sys- 
tems and for one dimensional systems ^0,01 (see 
also Q; reviews can be found in J, 10]). Semiclassical 
calculations involving complex trajectories in the mixed 
representation {x\e'''^^"'"/'^\z), on the other hand, were 
introduced in ll] and recently rediscovered 12] for the 
one-dimensional case (see [l^ for a different approach). 
Since the mixed representation is the most interesting 
for the propagation of wave packets, our purpose here is 
to generalize this formalism to many dimensions and to 
present some applications. 

The calculation of complex trajectories involves two 
difficulties: first, the effective dimensionality of the phase 
space is doubled, since both real and imaginary parts of 
position and momenta have to be computed; second, the 
trajectories must satisfy mixed boundary conditions, part 
at the initial time and part at the final time, a problem 
known as 'root search'. Therefore we also consider the 
possibility of employing only real trajectories in the semi- 
classical approximation. This is done by approximating 
the complex trajectories by real ones, satisfying modi- 
fied boundary conditions that are less restrictive than the 
original ones. Although such real trajectories approxima- 
tions are always less accurate than the original complex 
one, they are much simpler and sometimes have practi- 
cally the same accuracy Our approach, both with 
complex and real trajectories, does not involve integra- 
tions over initial conditions, a procedure that is common 
to IVR - Initial Value Representations (recent reviews 
of this method can be found in 14]). IVR methods are 
usually easy to apply and reasonably accurate for long 
times. Nevertheless, for short times the present method 
provides a much clearer physical picture since only a few 
families of trajectories are required. 



We start from a coherent state ]z), where 

z=i=(i?-lq-f*C-lp), (1) 

and the d-dimensional vectors q and p arc the average 
values of position and momentum for this state. The 
diagonal matrices B and C contain the position and mo- 
mentum uncertainties, respectively, and satisfy the con- 
dition B = hC~^. The position representation of this 
coherent state is a Gaussian, 

(xjz) - AAexp I ^p^(x - q) - i(x - q)^^-^ (x - q) | , 

(2) 

where Af — ]i?]~2 7r~4 (we use the symbol ] • ] for the 
determinant). After a time T, the propagated wave func- 
tion is given by 

^(x,r) = (xlif(T)lz), (3) 

where K{T) = e''"^/^'. 

In order to calculate the wave function semiclassically 
we shall follow the procedure of [ill IT^ . We first insert 
in 10) a resolution of unity to obtain 

^(x,r) = j dx'(xix(r)ix')(x'iz), (4) 

and substitute the quantity_6cL^(T)]x') by its semiclas- 
sical Van- Vleck expression [l^ [l^ . Then we make the in- 
tegration by the stationary exponent approximation. We 
shall see that the stationary points are in general com- 
plex numbers, and thus a deformation of the integration 
contour into the complex plane is unavoidable, taking the 
classical trajectories involved in the approximation to a 
complex phase space. 

The Van- Vleck formula in d dimensions is 

(xlA'(r)lx')yy = (27r*n)-'^/2 ^^^^gi^, (5) 

where ^(x, x',T) is the action of the classical trajectory 
that goes from x' to x in time T and S'xx' is the matrix 
of its second derivatives (we have incorporated Morse 
phases in S). If more than one such trajectory exists, 
one should sum their contributions. Before performing 
the integration, let us express the determinant in in 
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terms of the elements of the tangent matrix. As shown 
in the appendix, 



(x|if(T)|x'; 
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(6) 



Note that in the position representation nothing is said 
about the momentum of the corresponding classical tra- 
jectory, and therefore it is not necessary to introduce any 
complexification. In the coherent states representation 
the boundary conditions are too stringent as one tries 
to specify not only the initial and final points (and the 
time) but also the initial and final momenta 10] . 

In the next section we shall calculate the integral 
in the saddle point approximation, valid in the semiclas- 
sical limit. In section III we develop further approxi- 
mations that involve only real classical trajectories. We 
show some illustrative numerical applications in section 
IV and present our conclusions in section V. 



II. COMPLEX TRAJECTORIES 

In the semiclassical limit the wave function in Q) can 
be written as 

where 

^$ = + p^(x' - q)] - i(x' - q)^i?-2(x' - q). (8) 

We evaluate this integral by the usual saddle point 
method, which consists in expanding the exponent to 
second order around its stationary point Xq, while the 
prefactor is simply evaluated at this point. After per- 
forming the resulting Gaussian integration, this leads to 
the semiclassical approximation 
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where 
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and $0 = $(x,x[„T). 

The stationary point Xq, determined imposing the con- 
dition V'$|^, = 0, is given by 



B-\^^-ci) = iC-\p-Y>'^), 



(11) 



where Po(x, XQ,r) = — V'S*]^, . Both pp and Xg are in 
general complex numbers, and the whole classical trajec- 
tory therefore takes place in a complex phase space. It 
leaves Xq at time with the complex momentum Pq and 
arrives at the real position x at time T. The classical 



action S and the tangent matrix M will also be complex 
in general. 

Using the relation (see appendix) 



|<i>x'x' 



IBP 



(12) 



the final result may be written in terms of the tangent 
matrix as 



V'sc(x,r) 



:exp{-$o}. (13) 



This generalizes the one dimensional formula presented 
in [ill IT^ . to which it reduces for separable systems 
and that has proven to be accurate in the evolution of 
wave packets in many different systems. It is of course 
exact for the propagation of a d-dimensional coherent 
state in free space and in potentials up to quadratic (har- 
monic oscillator, charged particle in constant electromag- 
netic/grativational field). Differently from the Van Vleck 
approximation, the prefactor involves the square root of 
a complex number (remember that | • | is a determinant, 
not a modulus), and its phase must be determined dy- 
namically with the condition that for T = we have 
Mxx = 1 and Mxp = 0. 

The semiclassical approximation (|13ll depends on com- 
plex trajectories (q(i), p(i)) satisfying the boundary con- 
ditions 

-^(B-iq(0)+*C-ip(0)) =z, q(T)=x, (14) 

where we have used the fact that (|ll|l can be written 
B~ix^ + iC-ip^) = B-iq + iC-ip. The final value of 
the momentum is not restricted and will be complex in 
general. Following Klauder and Adachi 0, 0, we may 
write the initial condition as 



q(0)=q-|-a;, p(0) = p + iCfi-^a;, 



(15) 



where w = a + i/3 is a complex vector to be determined. 
The first condition is automatically satisfied for any w. 
For a fixed time T the propagation of this complex initial 
condition defines a complex map u) — s- q(T), the proper- 
ties of which have been studied in detail for the one- 
dimensional case in |^ . Only for some values of a; will 
it happen that q(T) £ W^, and we denote the set of all 
those points by fi. It is easy to see that w = belongs to 
il, in which case we have the classical trajectory of the 
center of the wave packet. 

However, the inverse of the map u) — > c\{T) is in gen- 
eral globally multivalued: there may be many trajecto- 
ries that end at the same q(T). Therefore Q, will con- 
sist in a finite collection of d-dimensional disjoint sets, 
called families. In the vicinity of a critical point (i.e., 
one for which dq{T)/dw is zero) the map is two-to-one, 
provided the second derivative is not zero. Such a crit- 
ical point is also called a phase space caustic. At these 



points |Mx 



-iMx 



0, thus preventing the validity of 
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the semiclassical calculation. It is possible to develop a 
semiclassical approximation based on the Airy function 
that remains valid near caustics. For the one-dimensional 
case this has been derived in [T^ . 

The family of trajectories that contains the point ui — 
is called the main family, and it provides the most im- 
portant contribution to the semiclassical approximation. 
As time increases, other families may become relevant. 
The imaginary part of $0 is positive for all trajectories 
that belong to the main family, but this may not be the 
case for other families. When Im($o) < one has a con- 
tribution that diverges when h ^ and therefore must 
be discarded. These are called non-contributing trajec- 
tories, and for some families it is necessary to introduce 
a cut-off in order to avoid them. 

Another delicate point is that of Stokes lines and expo- 
nential dominance, which is intrinsic to many asymptotic 
formulations. In the usual one-dimensional WKB, for 
example, the semiclassical approximation for stationary 
states becomes singular at classical turning points, and 
one must connect different local solutions by an analytic 
continuation. In so doing, one finds that there must be a 
change in the number of contributions along certain lines 
called Stokes lines [11 El IM EH ■ In the vicinity of such 
a line one contribution dominates exponentially over the 
other, and one is free to place a cut-off (the error due 
to the cut-off is less than the error due to the semiclas- 
sical approximation). The same phenomenon appears in 
the present formalism. Even though the location of these 
lines is hard to determine in principle, in practice when 
crossing them there appears a false divergence in the ap- 
proximation, which can be easily detected ^3- 



and each contribution to the propagator is obtained by 
a linearization of the dynamics. This method was ini- 
tially used to propagate wave packets and later to 



obtain coherent state correlation functions 



iHT/h\ 



in chaotic systems j2^,2Sj • The calculations we present in 
this section are close in spirit to these works, but instead 
of following the 'cellular' approach we start from the com- 
plex trajectory approximation and we also consider 
a variety of boundary conditions that the real trajectories 
may satisfy. Using different boundary conditions we ob- 
tain the 'central' trajectory approximation |23j | and also 
more general results similar to the 'off-centered' one pre- 
sented in psf. 

This section regards only calculation of wave functions, 
but a discussion of the quantity (z'|e~*^^/''|z) that pro- 
ceeds along the same lines may be found in |2q . 



A. Approximation via central trajectory 

The classical trajectory that starts at (q, p) will end, 
after a time T, in the point (qr,Pr)- Following ^3 we 
write 

x[, = q + Ax' (16) 
p'o = p + Ap' = p-S'x'x'Ax'-S'x'xAx (17) 
x = + Ax (18) 
p = p,. + Ap = pr + 5'xx'Ax' + SxxAx. (19) 

The stationary exponent condition can be written as 
Ap' — ihB~'^Ax', and equation lfT7|l can be solved to 
give (see appendix) 



III. APPROXIMATIONS BASED ON REAL 
TRAJECTORIES 

One may wish to find approximations for the expres- 
sion (|13|l that involve only real trajectories. There are 
many such possibilities. One possible choice is the tra- 
jectory that starts at the real point q with initial momen- 
tum Pi, different from p, and after a time T arrives at x. 
Another possibility is a trajectory that starts with mo- 
mentum p but from a different point q^ and also arrives 
at x. We can also give up the final point condition, for 
example, by choosing the unique trajectory that starts 
at q with momentum p. All these possibilities are sim- 
ilar to the ones already existent in the one-dimensional 
case [T^ . but in more dimensions one can in principle 
come up with others. For example, in two dimensions 
a trajectory may exist that starts at {qx,(lyi) with mo- 
mentum [pxiiPy) and ends at {x,y), but with qyi ^ qy 
and pxi 7^ Px- All these real trajectories should be good 
approximations for the complex stationary trajectory if 
the latter is not too deep into the complex plane. 

An important method that is also based on real tra- 
jectories is the 'cellular dynamics', developed by Heller 
|22l I23II , in which a grid of initial conditions is evolved 



Ax' = S(Mxx + iAfxp)"^S"iAx. (20) 

Now we expand the exponent in (|13|l around this tra- 
jectory to second order in Ax. The expansion of the 
action is 

S ~ Sr+p^Ax - p^Ax' 

+5<'^'=^-'>(5rJ;:0(A5)- 

The remaining terms are simply p-^Ax', which cancels 
out, and Ax'-^S^^Ax'. In the quadratic terms we intro- 
duce the tangent matrix and use (|20|l to obtain 

^'^P^^'^) - /L l-M | ^"P^ft^->' (22) 
where the exponent is given by (see appendix) 

1$, ^USr + P^Ax) - ^AxS-iSS-iAx, (23) 
h h 2 

where S = (-A/pp — iAfpx)(Afxx+*-^^xp) ""^^ Note that this 
is always Gaussian in Ax, with variable width. Therefore 
this approximation can never account for interferences 
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or tunnelling effects. Notice that while the formula H13() 
involves an infinite number of classical trajectories, at 
least one for each value of x, the one we just derived 
requires only the trajectory that starts in (q, p). For 
this reason this is called an Initial Value Representation 
(IVR). 

This formula was first derived by Heller [23 (see also 
psf) and is called the Thawed Gaussian Approximation 
or TGA (it was rederived with some detail in It 
becomes exact in the semiclassical limit h ^ (for a 
fixed value of time) and has been used, for example, in the 
study of decoherence 28] and of scars in quantum chaotic 
systems ^29j . In the applications presented here we are 
interested in quantum effects that cannot be reproduced 
by the TGA, and thus we do not consider it any further. 



B. Approximation via trajectory q ^ x 

Let us fix the initial coordinate of the trajectory, q, 
and demand that after a time T it arrives at x. We need 
to find the initial momentum pj for such trajectories, and 
in fact there may be more than one that satisfy the above 
conditions. We write 

x'o = q+Ax' (24) 
p'o = p, + Ap' = p, -S-x'x'Ax'. (25) 

Note that the complete expansion of Pq to first order 
should be p^ — 5x'x' Ax' — S'x'xAx, but we are keeping x 
fixed. 

Equation gives 

B-^Ax' ^iC-\Ap- Ap'), (26) 

where Ap = p — Pi . Using (|25|l we find 



The exponent contains the real action S'q and a term 
which is Gaussian in the difference between p^, the initial 
momentum of the trajectory, and p, the average momen- 
tum of the initial coherent state. It is important to note 
that Pi usually depends on x in a complicated manner, 
and thus the final wave packet will not, in general, be 
Gaussian. Also, there may exist more than one value 
of Pi, and a sum over all possible trajectories would be 
required, resulting in interference terms. Since the tra- 
jectory involved in the calculation depends on the initial 
q and final x points, this is not an IVR. 



C. Approximation via trajectory p ^ x 

We now fix the initial momentum of the trajectory and 
allow it to start from a point q^ that is different from the 
center of the wave packet. We write 



Po 



q, + Ax' (32) 
p + Ap'-p-5x'x'Ax', (33) 



and use the stationary exponent condition to find 



B-2)Ax' = -B-^Aq 



(34) 



or Ax' = — $^,^,_B^^ Aq, with Aq = q— q^. Once again, 
we expand the exponent in (|13() to second order in Ax', 
but this time we write it in terms of Aq. The final result 
is 



V'p(x,r) 



iM. 



xp 



exp{-$p}. 



(35) 



(■^5x'x' - B-2)Ax' = $x'x'Ax' = -^Ap, 



which we can invert to write 



A^' = -^*x'x'Ap. 



(27) 



(28) 



We now expand the exponent in H13() around this trajec- 
tory to second order in Ax'. Proceeding analogously to 
the one-dimensional case jl2j we obtain 

which, with a few algebraic manipulations, may be ex- 
pressed in terms of the tangent matrix (see appendix) as 



where 

l$p = l(5p-p-Aq) 

- iAq^B-i(Afxx + iM^p)-Hl^^B-'Aq. (36) 

We have obtained a Gaussian again, this time in the 
difference between q^, the initial position of the trajec- 
tory, and q, the average position of the initial coherent 
state. This is again not an IVR, and after a time T it 
will not result in a Gaussian in x. It may as well display 
interference between different existent classical trajecto- 
ries. 



= ~ ^Ap^C-i(Afxx + iMxp)-iMxpC-iAp. 

(30) 

The wave function becomes 

V'q(x,T) = .^^ , exp{-$q}. (31) 

^|Mxx + «A^xp| ^ 



D. Approximation via a mixed trajectory 

We now restrict ourselves to a two-dimensional system 
and choose the real trajectory that starts at [q^, qyj) with 
momentum {pxi,Py) and ends at {x,y), but with qyi ^ q-y 
and pxi 7^ Px- This time we have mixed conditions, and 
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we set 



q, + Ax\ (37) 
qy^ + Ay', (38) 

P'xO = Pxi + ^P'x = Pxi - Sx'x'Ax - Sx'y'Ay', (39) 



Xq 



PyQ 



Py 



Ap' ^Py- Sy'x'Ax' ~ Sy-y'Ay'. (40) 



Using these equations, the stationarity conditions can be 
cast in the form Ax' = — $~,^,A^, where 



AC - f i^Px/^ 



Ax' = 



(41) 



is S* « Sm + V'5''^Ax' + iAx''^S'. 
P^(x;,-q) = 



The expansion of the exponent to second order in Ax' 

Ax' for the action, 
V'S'-^Ax' — pyAy for the term involving 
the wave packet momentum and (xq — q)-'"i?^^(xQ — q) = 
Ax''^B-^Ax.'+b~'^Ay{Ay-2Ay') for the quadratic term. 
The hnear terms in Ax' cancel, and after we change from 
Ax' to A^ we have 



M 



-(SM-PyAy)--Ae^Z'.'^^ 



Ay2 
2bl 



(42) 



where Sm is the action for the mixed condition trajec- 
tories. Once again this can be written in terms of the 
tangent matrix, as we show in the appendix. The wave 
function becomes 



ipM{x,y,T) = 



exp{i$A/} 



(43) 



E. Alternative derivation 



Given the integral in |(7J) one may argue that, if the po- 
sition uncertainties are very small, only the region around 
q will be relevant for the integration. Expanding the ac- 
tion to second order around this point we have 



order in the difference x' — q (the same procedure was 
used in Isi^). After changing the integration variable in 
(I45|l from x to Pi, Vanicek and Heller arrive at a 
semiclassical result that is free of caustics. Even though 
expanding the action to first order only is inaccurate for 
simple systems such as the free particle and the harmonic 
oscillator, their final formula seems to work well in prac- 
tice. The expansion to second order we just presented 
is in principle more accurate, but the result it gives for 
the semiclassical fidelity is sensitive to caustics and thus 
probably less stable in numerical calculations. 

Finally we note that, if instead of inserting a position 
representation of unity in (jSjl, we used a momentum rep- 
resentation, 



V'(x,r) = J dp'(x|i^(T)|p')(p'|z) 



(46) 



then after a similar second order expansion of the action, 
this time around p' = p, we would arrive at the expres- 
sion H35|) for the wave function. This is justified when the 
momentum uncertainties are small. The TGA approxi- 
mation can also be obtained this way: one must enforce 
a stationary phase condition on the imaginary part of 
^(x, x') alone. 



IV. APPLICATIONS 

In this section we present a few numerical applications 
of the approximations we have just derived. We com- 
pare the complex-trajectories formula H13I) with the ones 
based on real trajectories, and also with exact quantum 
mechanical calculations, which we have carried out using 
Fast Fourier Transform methods. The purpose here is 
not to obtain extremely accurate numerical results, even 
though sometimes this is the case, but rather to illus- 
trate the usefulness of semiclassical calculations in many 
different situations. 



l<i>«-^5(x,q)--i(p,-p)^(x'-q) 
+ i(x'-q)^<i>x'x'(x'-q), 



(44) 



where p^ = — ^''5'lx'=a initial momentum, gener- 

ally different from p. Proceeding with the integration, 
we find the same result as in section V.A, which is based 
on the trajectory that starts at q with momentum p^ and 
ends at x with any momentum at time T. 

Jalabert and Pastawski |3QI| have used a similar argu- 
ment in their treatment of the quantum fidelity 



V/(x,T)^v(x,T)dx, 



(45) 



(in this equation ipi^j T) and V-'y(x, T) are obtained from 
an initial wave function by evolving it with two differ- 
ent Hamiltonians) , but they expanded the action to first 



A. Attractive Gaussian potential 

We start by investigating the semiclassical propagation 
in a two dimensional attractive Gaussian potential. 



V[r) 



(47) 



where = + y'^. A one dimensional version of this 
problem was already considered in [l^ . where the semi- 
classical approximation was shown to be very accurate. 
This potential is also interesting because, unless the par- 
ticle's momentum is very low (which is not the case we 
are interested in), there is only one classical trajectory 
that contributes to the real semiclassical formulas pre- 
sented in sections UFA to IIFD. In the complex case 
there may be more than one trajectory, but we will con- 
fine ourselves to the main family only, since it already 
gives a very good result. 
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We place the wave packet initially at q = (—10, 1), and 
chose bx = by = 1, so that the impact parameter is equal 
to the wave packet width. After a time interval of i = 4 
the main peak has followed a curved trajectory, arriv- 
ing at a negative value of y, and a smaller peak appears 
around j/ « 2, as we can see in Fig la. This is accurately 
reproduced by the semiclassical approximation ■(/'sc(x, t), 
shown in Fig. lb. The secondary peak is recovered al- 
most exactly, but the height of the main peak is wrong 
by a factor of 2 (notice the particular scale that has been 
used) . The phase of the wave function is also recovered, 
and in fact the overlap 

|(^|^.c)P (48) 

is around 92%. It is important to notice that the discrep- 
ancy comes from a small region around the peak, and 
that the functions agree very well at all other points. We 
have also calculated the real trajectory approximations 
■0q(x, t), V'p(x, t) and ■i/'M(x, i), but none of them can be 
distinguished from the complex one at this scale. 

The erroneous increase in the main peak is probably 
due to the presence of a caustic in complex phase space. 
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FIG. 1: (color online) Exact (up) and semiclassical (down) 
probability densities (times lO'^) at T = 4, with q = ( — 10, 1) 
and p = (3, 0), in the case of an attractive Gaussian potential. 
Except for the main peak (notice the change in scale), the 
wave function is accurately reproduced, and KV'IV'sc)!^ ~ 92% 



Even though only one real trajectory exists, in the com- 
plex case there may be more than one, leading to crit- 
ical points in the map described in section II. In order 
to obtain a better approximation in the vicinity of the 
peak, either a uniform approximation or incorporation 
of this secondary family of trajectories would be neces- 
sary. Finding these trajectories in practice is the noto- 
rious root-search problem, known to be very difficult in 
more than one dimension. The accuracy of the simpler 
formulas H31|l . (|35|) and H43|l in this case shows that they 
can be of practical use. 

B. A bound system 

We now study a bound system, an isotropic quartic 
oscillator: 

Vir) = Ar^ + Br'^, (49) 

where A — 0.5 and B — 0.1. The initial wave packet 
has parameters bx — by ~ 1, q — (0,0) and p = (2,0), 
which corresponds to a classical initial condition that is 
periodic with period r « 4.7. In Fig. 2a we show the 
probability density at T = 2.4, approximately half the 
classical period. It has a main peak at the origin and a 
small shoulder around x w —1.5. 

It is interesting that the approximation ipq{x, y, T) be- 
comes discontinuous in this case, as we see in Fig. 2b. 
This happens because only a certain region of coordi- 
nate space can be reached by real trajectories that start 
in the initial point q with an initial momentum that is 
close to p. Points outside of this region can eventually be 
reached, but the initial momentum must be so different 
from p that the actual contribution to ipq is negligible. 
In the border of this region there is a caustic, where the 
wave function diverges, and in the numerical simulations 
we must make this region a little smaller in order to avoid 
this. All approximations based on real trajectories suffer 
from this shortcoming, except for the IVR ipqp, which is 
always Gaussian. 

The approximation ipse, on the other hand, is based on 
complex trajectories and is well behaved in this case. It 
is presented in Fig. 2c, where we can see that it repro- 
duces very well the main peak. In fact, its only defect 
occurs near the shoulder. This is so because we have used 
only the main family, and in that region a contribution 
from a secondary family should be taken into account in 
order to give a good approximation (a similar effect can 
be observed in a one dimensional quartic oscillator [12| . 
where finding the secondary family is much easier). It 
is interesting that in this case ipse is far superior to the 
simpler real trajectory approximations, in contrast with 
the previous example where no caustics appeared, and 
the overlap H48|l in this case is around 95%. 

A better picture of the behavior of these wave func- 
tions is given in Fig. 3, where we show a cut along the 
line y = of the previous plots. The exact probability 
density is the solid line, while \tpsc(x,y)\'^ is the dashed 
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FIG. 2: (color online) Probability density at T = 2.4 for 
the quartic oscillator, with q = (0,0) and p = (2,0). The 
upper panel shows the exact calculation, the middle one is 
|i/)q(a;, i/)p and the lower one is \ipac{x,y)'^ . Using only real 
trajectories the result is very poor, but it becomes excellent 
when complex ones are used: the overlap 14811 in this case is 
around 95%. 



line and |-0q(a;,?/)p the dotted line. The first two agree 
well except around the region where the exact calcula- 
tion has a small shoulder. Inclusion of other families 
would certainly improve this result. As already noted, 
the approximation based on the real trajectory q ^ x 
fails completely for positive x because of the presence of 



0.30 




FIG. 3: Cut of the probability densities in Fig. 2 along the 
line y = G. The solid line is the exact result, the dashed line 
is |i/)ac(a:, t/)P and the dotted line is |i/)q(a;, j/)^. Notice that 
the latter must be cut because of the presence of a caustic. 



a caustic line. 



C. Circular billiard 

As our third example, we consider the motion inside 
a circular billiard with hard walls. If the particle is ini- 
tially at the center of the circle, the classical trajectories 
and also the tangent matrix can be computed analyti- 
cally, and we therefore consider this case only. An exact 
calculation for T = 0.5 is presented in Fig. 4a, where 
we have used p = (4, 0) and the radius of the billiard is 
i? = 3 (once again we use hx = hy = 1). As the packet 
approaches the wall, it develops interference fringes in 
the radial direction. 

We consider only the real approximation T/^q, but in 
this case for all final points x we should take into account 
the contribution of the many trajectories that reflect at 
the boundaries of the billiard. The actual number of such 
trajectories is infinite, but we consider only the two short- 
est ones, respectively with zero and one reflection, which 
give the main contributions. This gives origin to interfer- 
ence, as we can appreciate from Fig. 4b. The agreement 
with the exact result is excellent: the curvature is prac- 
tically the same, as well as the height and the position of 
the peaks. It is important to note that there is a collision 
with a hard wall involved, and thus an extra phase of 7r/2 
must be introduced in the contribution of the reflected 
trajectory. The overlap H48|l in this case is around 97%. 
Since this approximation is already very good, we do not 
present the complex calculation. We show again a cut 
along the line y = of the probability densities, in Fig. 
5. The small discrepancy could be corrected if a twice 
reflected trajectory was included. 
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D. Tunnelling system 

Finally we consider a system in which the tunnel effect 
plays an important role. We take a potential of the type 

V{x,y) = yoexp{- ^'''~;°^' }, (50) 

with = + j/'^, which describes a circular ridge of 
radius ro in the plane, centered around the origin. When 
an incident wave packet with energy less than Vq is scat- 
tered by this potential, there is a probability that the 
particle will tunnel into the ridge. In Fig. 6a we see 
the exact calculation ai T — 2.5, for a potential with 
Vb = 10, ro = 5,(T = 10 and an initial wave packet with 
q = (-10,0) and p = (4,0). The total probability of 
being located inside the ridge is around 10% in this case. 

A semiclassical calculation for tunnelling though a 
square barrier involving complex trajectories was pre- 
sented in , where only the coherent state representation 
{z'\e~^-^^'^\z) was considered. In the present case the 
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FIG. 4: (color online) Probability density at T — 0.5 in the 
case of a circular billiard, with q = (0, 0) and p — (4, 0). The 
upper panel shows the exact calculation and the lower one is 
|'!/'q(a;, J/) p. Using only real trajectories we have a very good 
result (KV'IV'q)!^ ~ 97%), including effects due to curvature 
and interference. 



0.5 




X 

FIG. 5: Cut of the probability densities in Fig. 4 along the 
line y = 0, displaying the exact (solid) and the semiclassical 
(dashed) results. The latter is obtained from the interference 
of a direct and a reflected trajectory. 

classical motion must be solved numerically and the pres- 
ence of turning points leads to the appearance of caustics. 
Nevertheless, provided the probability amplitude is not 
large in the vicinity of the caustics, the real trajectory 
approximation |-0q(a;, y)^ is able to give an accurate re- 
sult, as we can see in Fig. 6b (the overlap between the 
transmitted wave function in the exact and semiclassical 
calculations is around 94%). This is easy to understand 
if we remember that for each value of the pair {x, y) we 
need a different initial momentum pi and, even though 
a classical particle with the average momentum p would 
be reflected by the potential, there will be values of pi 
for which transmission is possible. The other real trajec- 
tory approximation |'0p(2;, y)Pi on the other hand, works 
poorly in this case, because it involves variation only on 
the initial position and this does not affect the energy of 
the trajectories. 

The full complex semiclassical calculation would give 
even better results than Fig. 6b, but this requires extend- 
ing the potential to the complex plane. This extension 
involves trigonometric functions that make the numeri- 
cal evolution very demanding. It is clear that the sim- 
plicity of the trajectories involved in the calculation of 
|V'q(a;, y)\^ is of great practical advantage. 

V. CONCLUSION 

We have generalized the semiclassical approximation 
for the propagation of wave packets based on complex 
trajectories derived in 0, 0| to multidimensional sys- 
tems. Several further approximations based on real tra- 
jectories were also derived from this basic formula, in par- 
ticular Heller's Thawed Gaussian Approximation (TGA). 
Apart from the TGA, all other formulas are not Initial 
Value Representations and are able to accurately repro- 
duce non-Gaussian wave functions and also quantum in- 
terference when more than one family of trajectories is 
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tion times. For long times the number of trajectories in 
bound systems increases and caustics proliferate, making 
a practical application of the formulas more difficult. If it 
is possible to overcome this problem, the study of chaotic 
systems would naturally be the next step. 
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FIG. 6: Probability density at T = 2.5 for the ridge potential, 
with q = ( — 10,0) and p — (4,0). The upper panel shows 
the exact calculation and the lower one is |i/)q(a;, j/)^ (times 
10^). Using only real trajectories it is possible to accurately 
reproduce tunnelling effects (|(i/'l''/'q)P ~ 94%). 



present. 

These theoretical results were tested in very distinct 
particular cases, starting with scattering by an attrac- 
tive potential, where the classical trajectories must be 
computed numerically. For positive energies this poten- 
tial has no turning points and thus no caustics. The 
complex and real approximations give indistinguishable 
results that are very close to the exact calculation. The 
second case was a bound nonlinear system, where a large 
number of contributing classical trajectories exist. Us- 
ing only the main family we obtained a very good result 
with the complex approximation. In this case the real 
trajectories approximations are not practical because of 
the many caustics involved. We also studied the mo- 
tion inside a circular billiard, taking into account two 
real trajectories for ipq(x,T), which displayed effects of 
curvature and interference. Finally, we considered the 
tunnel effect and showed that again '0q(x, T) is able to 
accurately reproduce the quantum result. 

All cases studied in this paper are integrable and have 
circular symmetry, which clearly introduces simplifica- 
tions. We have also considered relatively short propaga- 



APPENDIX 

Consider a classical trajectory, satisfying Hamilton's 
equation 



di 



JWH, 



(A.l) 



where J is the usual symplectic matrix and V is the 2c?- 
dimensional gradient. A variation around this trajectory 
satisfies 



d_ 

dt 



J 



pp 



(A.2) 



where the second derivatives of H are computed at the 
reference trajectory. Multiplying both sides on the left by 
a matrix containing the inverse quantum uncertainties, 
B and C, and inserting an identity in the r.h.s. we can 
rewrite (IA.2|I as 



d_ /<5x 

dt Up 



B ^HpxB B ^HppC 



6i 
Sp 



(A.3) 



where x = B^^x and p = C^^p. 

Now consider a trajectory that starts from x' with mo- 
mentum p' and arrives at x with momentum p (not re- 
lated with the initial coherent state label), and suppose 
we make small displacements in its initial and final co- 
ordinates. This induces variations in the initial and final 
momenta according to 



Sp 
Sp' 



>5'xx 
^•^x's 



^x 



5x 



(A.4) 



On the other hand, the tangent matrix is defined to be 
the linear application that relates the initial and final 
displacements. 



S5i 
Sp 



Mxx Mxp 



dx' 
Sp' 



(A.5) 



where we have included explicitly the quantum uncer- 
tainties for convenience. Inverting equation (jA.4|l it is 
possible to show that 



^x'x|-p|A^xp| 



(A.6) 
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which we have used in equation l^. It is also possible to 
show that 5x'x' — CM^^ M^^ixB^^ and therefore 

$x'x' = B-l(^^4-plMxx - l)B-\ (A.7) 

where we have used C/h = B^^ . This leads to 

i"^ |Mxx + iMxp| 



l^x'x'l 



ISP lAf. 



(A., 



as stated in equation The inverse of <i>x'x', used 

in H43() . can also be expressed in terms of the tangent 
matrix: 



^•x'x' - -^B{M^^ + iMxp)-'M^pB. (A.9) 



If we now take the time derivative of equation (jA.Sfl . 
and compare the result with HA.3() we conclude that 



dM 
dt 



-C^^Hy-y-B 



B ^HppC 



-c- 



M. 



(A.IO) 



This is the dynamical equation for the tangent matrix, 
which may be simplified for the large number of cases 
in which iJxp = ^^px ~ and Hpp is the inverse of 
the mass. In practical applications these may be solved 
together with the equations of motion, making it possible 
to follow the phase of the prefactor in H13(l . 
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